Supersolid polar molecules beyond pairwise interactions 
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We explore the phase diagram of ultracold bosonic polar molecules confined to a planar opti- 
cal lattice of triangular geometry. External static electric and microwave fields can be employed 
to tune the effective interactions between the polar molecules into a regime of extended two- and 
three-body repulsions of comparable strength, leading to a rich quantum phase diagram. In addition 
to various solid phases, an extended supersolid phase is found to persist deep into the three-body 
dominated regime. While three-body interactions break particle-hole symmetry explicitly, a charac- 
teristic supersolid-supersolid quantum phase transition is observed, which indicates the restoration 
of particle-hole symmetry at half-filling. We revisit the spatial structure of the supersolid at this 
filling, regarding the existence of a further supersolid phase with three inequivalent sublattices, and 
provide evidence that this state is excluded also at finite temperatures. 

PACS numbers: 67.80.kb, 75.40.Cx, 64.70.Tg, 75.40.Mg 



I. INTRODUCTION 

Identifying conditions for bulk supersolidity, both in 
the continuum and for bosons confined in periodic po- 
tentials, is a major task in order to explore the super- 
solid state of matter i — A fascinating approach towards 
probing the range of stability and the nature of super- 
solids would be provided by the direct manipulation of 
the inter-particle interactions. In this respect, ultracold 
polar molecules^ — exhibit the remarkable opportunity to 
tailor many-body interactions by the application of static 
electric and microwave fields. 8 10 It has been demon- 
strated, that by appropriately tuning the external control 
parameters, polar molecules can in fact be driven into 
a regime, where three-body repulsions strongly compete 
with effective two-body interactions between the parti- 
cles^ Starting from a supersolid phase, this tune-ability 
thus allows to explore the stability of lattice supersolid 
states beyond the regime of pairwise (two-body) interac- 
tions. 

Here, we consider in particular a two-dimensional 
setup, where the molecular dipole moments align par- 
allel to a static electric field directed perpendicular to 
the confining plane. The effective dynamics of the polar 
molecules in an optical lattice potential is then described 
by a Bose-Hubbard Hamiltonian of hard-core bosons 
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where b{ and b\ denote boson annihilation and creation 
operators respectively on lattice site z, and the local 
density operator rii = b\bi. Furthermore, t denotes 
the nearest-neighbor hopping matrix element, and \i the 
chemical potential, controlling the filling n of the lattice. 
In contrast to two-body terms, the three-body interac- 
tions explicitly break the particle-hole symmetry that the 
model exhibits in the limit W = 0. The interactions 
furthermore depend on the lattice geometry via their re- 
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FIG. 1. (Color online) Leading contributions to the many- 
body interactions for polar molecules on a triangular optical 
lattice. Numbers below the interaction terms denote the rel- 
ative strength of the two- (three-) body terms in units of the 
nearest-neighbor values V (W). 



spective strengths, 
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where denotes the distance between lattice sites i and 
j. We consider here the case of a triangular lattice, for 
which the geometric factors of the leading contributions 
to the interaction terms are given in Fig. 1. Analyzing 
the dependence of these parameters on the experimental 
setup shows that W/V can be tuned up to about 0.7^ 
Furthermore, t can be varied independently by control- 
ling the depth of the optical lattice. 

We explore the consequences of three-body interac- 
tions on the low-energy properties. Of particular inter- 
est are the stability and nature of the supersolid (SS) 
state of bosons on the triangular lattice found in the 
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FIG. 2. (Color online) Ground state phase diagram of polar 
molecules on a triangular lattices. The upper (lower) panel 
shows the regime W > V (V > W), at t/W = 0.1 (t/V = 
0.1). SF denotes superfluid regions, SS-A and SS-B the two 
supersolid regimes, that meet at a first-order quantum phase 
transition along the line of half- filling, n — 1/2 (squares). 
Diamonds trace the SS-SF transition line. The inset focuses 
in on the large- W region. 



particular cases of nearest-neighbor— ~— as well as for 
1/rfj -extended two-body repulsions. 17 For the nearest- 
neighbor case, a discontinuous quantum phase transition 
separates two distinct SS phases, which differ in their in- 
ternal structured In the SS-A (SS-B) phase, of density 
n above (below) 1/2, two sublattices of the triangular 
lattice show density fluctuations about a larger (smaller) 
common mean value than the sites of the third sublat- 
tice. These phases are related by the exchange of parti- 
cles and holes, and meet at the particle- hole symmetric 
line of half-filling, where another possible SS state (SS- 
C), with distinct occupations on all three sublattices 
is indeed excluded in the ground state by the global U(l) 
symmetry of particle number conservation^ 

In order to study the effects of three-body repulsions 
on the SS phases, we first restrict to the nearest-neighbor 
two- and three-body terms in H, i.e. the most basic 
model with competing two- and three-body interactions 
on the triangular lattice. Then, we assess the effects of 
longer-ranged interactions. The stability of the SS phase 
under three-body interactions is observed already within 
the reduced Hamiltonian. This is unlike the case of the 
square lattice^ for which the emergence of a SS state 
relies on interaction ranges beyond nearest-neighbors. 

II. METHOD 

We employed quantum Monte Carlo simulations us- 
ing a generalized directed-loop algorithm in the stochas- 
tic series expansion representation* 2 ^— Decoupling the 
Hamiltonian in terms of 6-sites clusters as in Fig. 1, 
allows us to account efficiently for all interactions of 



strengths larger than 0.016F, and 0.05W, respectively. 
The simulations were performed on finite systems with 
periodic boundary condition, containing N = 1? lattice 
sites (L denotes the linear system size), with L up to 
48. In order to obtain ground state properties for these 
finite systems, the temperature T needs to be chosen suf- 
ficiently low, as specified below. 



III. PHASE DIAGRAM 

The ground state phase diagram for the nearest- 
neighbor model is shown in Fig. 2 in the vicinity 
of the supersolid phase near half-filling and for t = 
0.1 max(VK, V). Besides the superfluid (SF) regime two 
solid phases appear of densities n = 1/3 and 2/3, respec- 
tively. The corresponding density structure factor 

relates to the density correlation function at the cor- 
ner of the hexagonal Brillouin zone at momentum Q = 
(27r/3,0)^ In addition, the phase diagram exhibits an 
extended SS region, where both S/N and the SF den- 
sity ps scale to finite values in the thermodynamic limit. 
Here, 

P5 = |(^ 2 ) (5) 

is obtained from measuring the boson winding number 
W fluctuations* 2 ^ The phase boundaries in Fig. 2 result 
from data such as those shown in Fig. 3, after an extrap- 
olation to the thermodynamic limit (cf. the inset of Fig. 
3 for an example). 

We find, that the SS state, realized on the triangular 
lattice in the absence of three-body interactions, remains 
stable even for large values of W/V . Only in the limit 
of purely three-body interactions (V = 0) is the SS de- 
stroyed. The solid phase of filling n = 1/3 is stable up to 
W ~ 2V. This destabilization of the n = 1/3 solid phase 
via W is consistent with results obtained for the classical 
limit (t = 0) of our model, i.e. the Ising model with two- 
and three-body interactions in a magnetic fiek L 24 i 25 



IV. EXTENDED MODEL 

We now turn to the extended model approximating 
the long-range nature of the dressed dipolar interactions 
by taking into account all interaction vertices drawn in 
Fig. 1. First, we consider the parameter region, where 
the SS phase is stabilized. In particular, Fig. 4 shows 
n, ps and S as functions of p at t/V = 0.1 for different 
ratios W/V < 1 around the SS regime. An extended SS 
phase is clearly identified between the n = 1/3 and the 
n = 2/3 solid. 
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FIG. 3. (Color online) Left: Filling n, superfluid density ps, 
and density structure factor S as a function of the chemical 
potential p for various values of W/V at t/V = 0.1, taken for 
L = 12 and T/t = 0.1. Right: Filling n dependence of the 
superfluid density ps at W/V = 0.3 and t/V = 0.1, taken for 
L = 12 (squares) and 24 (solid circles) and T/t = 0.1. Inset: 
Finite-size extrapolation of ps (circles) and S/N (squares) at 
n = 1/2. 
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FIG. 5. (Color online) Density n as a function of the chemical 
potential p within the extended model for different values of 
t/W at V = 0, taken for L = 12 and T/t = 0.1. 



these solids are shown in Fig. 6. 
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FIG. 4. (Color online) Filling n, Superfluid density ps, and 
the density structure factor S as functions of the chemical 
potential p for different values of W/V at t/V = 0.1, taken 
for L = 12 and T/t = 0.1 in the extended model. 



Next, we explore the phase diagram of the extended 
model beyond the SS region. Of particular interest are 
incompressible phases that are stabilized by the three- 
body interactions. To explore these, we first consider the 
limiting case of purely three-body interactions (V = 0). 
Fig. 5 shows quantum Monte Carlo data of the filling as 
a function of p for different ratios t/W at V = 0. In 
addition to the previously identified density plateaus at 
fillings n = 1/3 and n = 2/3, incompressible phases at 
fillings n = 3/8, 1/2, and 3/4 are seen to emerge upon 
decreasing t/W. The spatial structures that characterize 



In the n = 3/8 solid, the system forms irregular pat- 
terns based on nearest neighbor dimers. From inspecting 
various real space density distributions from the quan- 
tum Monte Carlo configurations, we conclude that these 
dimers are not arrange in any periodic pattern. This 
dimerization leads to a strong energy penalty for finite 
values of V, such that this phase is unstable for V > W. 
Similarly, the n = 1/2 phase, where stripes of empty and 
occupied lattice sites alternate, is not observed within 
this parameter regime. The solid phase at n = 3/4 
however remains stable; for example we observe it for 
W/V = 0.2, t/V = 0.1 around p/V = 14 (not shown). 

Due to the frustration induced be the competition be- 
tween the two- and three-body interaction terms, we 
are able to perform ergodic simulations for the extended 
model only for ratios W/V up to about 0.3. For example, 
as the dip in the data for S/N near p/V = 5.5 in Fig. 4 
illustrates, we cannot reliable resolve the transition re- 
gion between the SS phase and the n = 2/3 solid for the 
extended model. Using annealing techniques, the simu- 
lations did not improve within the transition region. It 
remains possible, that another, maybe incommensurate, 
phase appears here. However, from our simulations, we 
are not able to conclude on this issue. Multi-body in- 
teractions at larger distances, as seen above for the case 
of purely three-body interactions, are in fact expected to 
stabilize additional plateaus at further fractional fillings, 
with the possible proliferation of low-energy metastable 
states^ and Spivak-Kivelson bubble-type transitions be- 
tween solid and superfluid regions^ However, a reliable 
numerical study of these issues is beyond the present ap- 
proach^ 




FIG. 6. (Color online) Real-space structures of incompress- 
ible phases stabilized in the extended model, labeled by their 
filling fraction n. For crystalline states, the crystal unit cell is 
also indicated. For the n — 3/8 phase, which does not exhibit 
a periodic structure, a typical state is shown, consisting of an 
irregular arrangement of isolated dimers (pairs of neighboring 
occupied sites). 
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FIG. 7. (Color online) Histogram P of n + for different sys- 
tem sizes taken at t/V = 0.1, W/V = 0.3, and T/t = 0.05 
within the SS-A and SS-B phase (left panel) and at the SS-SS 
transition at n = 1/2 (right panel). 



V. SS-SS TRANSITION 

In the curves of the superfluid density as functions of 
ji in Fig. 3 (left panel), dips occur for values of p/V near 
3.1. This suppression of ps appears around n = 1/2, 
as seen from Fig. 3 (right panel), which shows ps as a 
function of the filling n. As seen from the finite size ex- 
trapolation in the inset of Fig. 3, the SS however remains 
stable in the thermodynamic limit at this filling. The 
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FIG. 8. (Color online) Histograms P of n + for different 
system sizes taken for the extended model at W/V — 0.1, 
t/V = 0.1 and T/t = 0.1 at half-filling n = 1/2. 



suppression of ps relates to a remarkable feature of su- 
persolidity of hard-core bosons on the triangular lattice: 
the presence of a SS-A to SS-B quantum phase transi- 
tion at half- filling. 16 A robust means of identifying this 
SS-SS transition is to monitor the histogram of the local 
densities.— Denoting by 
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(6) 



the imaginary-time averaged local density at lattice site 
i for a given Monte Carlo configuration, one obtains the 
fraction of sites with mean occupation fi{ > 1/2 from the 
observable 
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where denotes the Heaviside function. Supersolid 
ground states belonging to the SS-A (SS-B) type are 
characterized by peaks in a histogram P of n + near 2/3 
(1/3), such as those in the left panel of Fig. 7. The ob- 
servation of a doubly-peaked histogram in the quantum 
Monte Carlo data exhibits the presence of both SS-A as 
well as SS-B configurations at n = 1/2, whereas the SS- 
C phase would correspond to detecting a singly-peaked 
histogram at n + = 1/2. Performing such an identifica- 
tion of the nature of the SS state for W ^ 0, we find 
that the suppression in ps is due to the SS-SS quantum 
phase transition at n = 1/2. As an example, we show 
in the right panel of Fig. 8 histograms taken for different 
system sizes at n = 1/2 for W/V = 0.3 and t/V = 0.1. 
One clearly identifies for sufficiently large system sizes a 
characteristic two-peak structure, indicative of the SS-A, 
SS-B coexistence at half-filling. This feature is not re- 
stricted to the nearest-neighbor model. In Fig. 4, for the 
extended model, similar dips are observed in ps as for the 
nearest-neighbor case. The histogram of n + , measured 
for the extended model at W/V = 0.1 and t/V = 0.05 
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FIG. 9. (Color online) Superfluid density ps (upper panel) 
and structure factor S (lower panel) at V/W = 0.5, t/W = 0.1 
for different system sizes. The different phases are highlighted 
(light) yellow for (the low density solid) SS-A and (light) 
red for (the high density solid) SS-B, respectively. The in- 
set shows finite size extrapolations for the structure factor at 
fi/W = 1.09 (a), 1.148(b), 1.438 (c), and 1.496 (d). 



and shown in Fig. 8, indeed also exhibits the characteris- 
tic double-peak structure of the SS-A, SS-B coexistence. 
Hence, although particle- hole symmetry is explicitly bro- 
ken for W 7^ 0, the characteristic SS-SS transition indi- 
cates its dynamical restoration in the SS region along the 
line of half- filling, n = 1/2. 

While the n = 1/3 solid is not stable for large values 
of W/V > 2 (cf. the inset of Fig. 2), we observe that 
within the remaining SS stripe both the SS-A and the 
SS-B phases remain stable. In fact, upon increasing p at 
V/W = 0.3, one enters the SS region from the SF region 
at n = 0.48(1) < 1/2, i.e. into the SS-A phase, before 
the transition to the SS-B phase takes place at n = 1/2. 
This stability of the SS-A phase furthermore leads to a 
SS re-entrance phenomenon when the density increases 
at a fixed ratio of V/W = 0.5: Fig. 9 shows the superfluid 
density p s and the structure factor S/N at V/W = 0.5 
and t/W = 0.1. Upon increasing the chemical potential 
/i/W, the system's ground state changes (i) from the 1/3- 
solid to a SS-A state at p/W w 0.9, then (ii) becomes 
a uniform superfluid at p/W w 1.2, before (iii) a transi- 
tion occurs around p/W ~ 1.5 to the SS-A phase, which 
eventually (iv) at fi/W ~ 1.9 enters into the 2/3-solid 
phase. In particular, we thus observe here re-entrance 
behavior with respect to the SS-A phase upon increasing 
the filling. This is indicated by the dashed line in Fig. 2. 



VI. ASSESSING THE SS-C PHASE 

A recent self-consistent cluster mean-field approach re- 
ports that the SS-C phase can be stabilized by thermal 
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FIG. 10. Histograms P of n + taken at W = 0, t/V = 0.1, 
fjb/V = 3 for different temperatures and system sizes. 



fluctuation J£ At n = 1/2, the SS-C state was obtained 
at t/V = 0.1 for T/t > 0.5, while the SS melting tem- 
perature was estimated near T/t = 6. This tempera- 
ture regime is not consistent with the results in Ref. 16, 
where it was found, that the SS melts already around 
T/t = 0.5. We thus performed finite- temperature simu- 
lations in order to assess the proposed SS-C scenario. In 
our histograms (cf. Fig. 10) at W = and t/V = 0.1 a 
double-peak structure and a suppression at n + = 1/2 
emerged for sufficiently large system sizes at least up 
to T/t w 0.25, indicative of the SS-A and SS-B coex- 
istence at n = 1/2. The peaks broaden and shift towards 
n + = 1/2 upon increasing T. For T/t beyond 0.25, we 
cannot resolve a double-peak structure, instead a single 
broad peak around n + = 1/2 persists up to the largest 
system size considered. Since we find that thermal dis- 
order leads to a substantial weight near n + = 1/2, this 
observation does however not provide evidence for a SS- 
C state, but instead suggests that the double-peak struc- 
ture is beyond our numerical resolution for T/t beyond 
0.25. Within quantum Monte Carlo, the SS-C thus re- 
mains elusive. We also did not observe it for W ^ 0. 



VII. CONCLUSIONS 

We showed that polarized polar molecules on the tri- 
angular lattice exhibit a robust supersolid phase in the 
presence of strong three-body interactions. The broken 
particle- hole symmetry is recovered along the line of half- 
filling, as indicated by the persistence of a first-order 
supersolid-supersolid transition. Furthermore, based on 
our findings, we exclude the existence of a SS-C phased 
also at finite temperatures. Experimentally, the differ- 
ence between a normal liquid, commensurate solids and 
supersolid phases can be detected via time-of-flight im- 
ages. Density correlations can also be probed using Bragg 
spectroscopy^ Furthermore, detecting the spatial struc- 
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ture of the solid and supersolid phases is achieved with 
a local addressability of individual lattice sites; the feasi- 
bility of this method for atoms has indeed been demon- 
strated recently^ — For the future, we plan to analyze in 
detail the quantum phase transition between the super- 
solid and the superfluid phase, for which a recent cluster 
mean- field theory^ predicts that for the nearest-neighbor 
model this transition is first-order away from the particle- 
hole symmetric point at fi/V = 3, where instead a con- 
tinuous transition was obtained, consistent with previous 
findings j 13 i 15 It will be interesting to verify this prediction 
and to clarify the unversality class of the later continuous 



transition. 
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